/*
  File:    TumorGrowthInitialization.cpp
  Purpose: Routines to initialize the tumor growth data. Includes 
           auto-detection of resume data. 

  Contact: Paul Macklin
           pmacklin@math.uci.edu
		  http://math.uci.edu/~pmacklin
*/

#include <cstdio>
#include <cstring>
#include <cctype>
#include <iostream>
#include <cmath>
#include "../Edinburgh_Angio/myangio.h"

using namespace std;

#include "../Matrix.h"
#include "../FileNames.h"
#include "../ParseText.h"
#include "../ImprovedMath.h"
#include "../BMPtoMatrix.h"
#include "../LogfileEntry.h"

#include "../Plotting.h"

#include "GaussianFilter.h"
#include "DiscreteReinit.h"
#include "Oxygen.h"
#include "Pressure.h"
#include "reinit.h"
#include "TumorGrowthInitialization.h"
#include "Curvature.h"
#include "TumorSaveUtilities.h"
#include "TAF.h"
#include "Pressure.h"
#include "Microenvironment.h"

#include "Polynomial.h"

#include "../tissue_types.h"

// This overwrites the target image with black where the 
// two images, and white otherwise. 

bool ImageDifference( BMP& Reference, BMP& Target )
{
 RGBApixel tempWHITE; 
 tempWHITE.Red = 255;
 tempWHITE.Green = 255;
 tempWHITE.Blue = 255;
 
 RGBApixel tempBLACK;
 tempBLACK.Red = 0;
 tempBLACK.Green = 0;
 tempBLACK.Blue = 0;

 for( int j=0 ; j < Reference.TellHeight() ; j++ )
 {
  for( int i=0 ; i < Reference.TellWidth() ; i++ )
  {
   if( Reference(i,j)->Red != Target(i,j)->Red ||
       Reference(i,j)->Green != Target(i,j)->Green ||
       Reference(i,j)->Blue != Target(i,j)->Blue )
   { *Target(i,j) = tempBLACK; }
   else
   { *Target(i,j) = tempWHITE; } 
  }
 }

 return true;
}


bool InitializeTumorGrowth( void )
{
 extern double _EPS;

 extern bool ResumeFromPreviousData;
 extern bool InitializeSimulationFromPictures;
 extern bool ScratchSpaceMode;
 
 extern double OutputInterval;
 extern double NextOutputTime; 
 
 extern double SimulationStartTime;
 extern double SimulationEndTime;
 extern double OutputInterval;
 
 extern double TumorOxygenUptakeRate;
 extern double QuiescentOxygenUptakeRate;
 extern double NecroticOxygenUptakeRate;
 extern double ECMoxygenUptakeRate;
 extern double ExteriorOxygenDiffusivity;

 extern char SaveDirectory [1025];
 extern char* ParameterFile;
 extern char Xfile [1025];
 extern char Yfile [1025];
 extern char PhiFile [1025];
 extern char PhiNecroticFile [1025];
 extern char PhiQuiescentFile [1025];
 extern char OxygenFile [1025];
 extern char OxygenDiffusivityFile [1025];
 extern char PressureFile [1025];
 extern char MobilityFile [1025];
 extern char ECM_File [1025];
 extern char MDEfile [1025];
 extern char TAFfile [1025];
 extern char TumorParameterFile [1025];
 extern char BloodVesselIndicatorFile [1025];
 extern char BloodVesselPressureXFile [1025];
 extern char BloodVesselPressureYFile [1025];
 extern char InterPressureFile [1025];
 extern char CellColorFile [1025];
 extern char LymIndicatorFile [1025];
 extern char PermeabilityFile [1025];
 extern char WallGradientFile [1025];
 extern char ConductFile [1025];
 extern char SproutTipsFile [1025];
 extern char StructureFile [1025];
 extern char BloodVesselDensityFile [1025];

 extern char TumorPictureFile [1025];
 extern char TissuePictureFile [1025];
 
 extern char TissueTypesFile [1025];
 
 extern char HaematocritLevelFile [1025];
 extern char BloodVesselIndicatorXFile [1025];
 extern char BloodVesselIndicatorYFile [1025];
 extern char BloodVesselIndicatorFile [1025];
 extern char SproutAgeFile[1025];
 extern char VesselAgeFile[1025];
 extern char BloodVesselRadiusXdirFile [1025];
 extern char BloodVesselRadiusYdirFile [1025];
 extern char BloodVesselRadiusFile [1025];
 extern char ConductXFile [1025], ConductYFile [1025]; //conductivity

 
 extern double MaxMobility;

 extern bool TopBloodVessel;
 extern bool BottomBloodVessel;
 extern bool IsFirstOxygenUpdate;

 extern double bulk_blood_tissue_transfer_rate;
 extern double bulk_blood_transfer_cutoff_distance; 
 extern double QuiescenceOxygenLevel;
 extern double NecrosisOxygenLevel;
 
 extern double ApoptosisRate;        // A in all papers
 extern double TumorAggressivenessG; // G in all papers
 extern double TumorNecrosisRateGN;  // GN in all papers 
 extern double DRG_Effect;
 
 extern double TumorProliferationConstant; // b in the paper
 extern double TumorSurfaceTensionConstant;
 extern double ExteriorMobility;
 extern double Custom_Cutoff_Max; 
 extern double collapseConstant; 
 extern double TAFdiffusionConstant;
 extern double TAFprodctionConstant;
 extern double TAFdecayConstant;
 extern double EndothelialTAFuptakeRate;
 
 extern double MDEdiffusionConstant;
 extern double MaxMDEdiffusionConstant;
 
 extern double MDEprodctionConstant;
 extern double MDEdecayConstant;
 extern double MDE_DegradationConstant;
 
 extern double DRGprodctionConstant;
 extern double DRGDecayVessel;
 extern double DRGdecayConstant;
 extern double DRG_DegradationConstant;
 extern double D_DRG;
 extern double AdvectionConstant;
 
 extern double ECM_DegradationConstant;
 extern double ECM_ProductionConstant;

 extern double Xleft;
 extern double Xright;
 extern double Ybottom;
 extern double Ytop;
 extern int Xnodes;
 extern int Ynodes;
 extern Matrix X;
 extern Matrix Y;
 extern double dx;
 extern double dy;
 extern double dt;
 extern double t;
 
 extern Matrix Phi;
 extern Matrix PhiNecrotic;
 extern Matrix PhiQuiescent;
 extern Matrix OutPressure;
 extern Matrix Oxygen;
 extern Matrix Pressure;
 extern Matrix MDE;
 extern Matrix TAF;
 extern Matrix ECM;
 extern Matrix Original_ECM;
 extern Matrix Mobility;
 extern Matrix OxygenDiffusivity;
 extern Matrix MDEdiffusivity;
 extern Matrix Structure;
 extern Matrix BloodVesselDensity;
 extern Matrix Shuffle;
 
 extern Matrix BloodNodePressure;
 extern Matrix FlowX, FlowY;
 extern Matrix WssX, WssY, Wss;
	extern Matrix SmX,SmY,SwssX,SwssY,SpX,SpY,SsX,SsY,SvarX,SvarY,Sm,Swss,Sp,Ss,Svar;
 extern Matrix Stability;
 extern Matrix TracerX, TracerY, Tracer, TracerS, TracerT, TracerD, Debuger, ShearRate,SrtX,SrtY, Macrof, Leak;
 extern Matrix VisX, VisY; //viscosity
 extern Matrix ConductX, ConductY; //conductivity
 extern Matrix BloodVesselRadiusXdir, BloodVesselRadiusYdir;
 extern Matrix BloodVesselRadius;
 extern Matrix SproutTips;
 extern Matrix SproutAge;
 extern Matrix VesselAge;
 extern Matrix HaematocritLevel;
 extern Matrix BloodVesselIndicatorX;
 extern Matrix BloodVesselIndicatorY;
 extern Matrix BloodVesselIndicator;
 extern Matrix VEGFR;
 extern Matrix Notch;
 extern Matrix Delta4;
 extern Matrix VesselCells;
 extern Matrix PressureGradientX;
 extern Matrix PressureGradientY;
 extern Matrix CollapseClock;
 extern Matrix HomPressure;
 
	extern Matrix TotalPressure;
	extern double VFWater;
	extern double VFCell;
	extern double ChMobility;
	extern double ChLength;
	extern double MitosisRate;
	extern double waterpressureA;
	waterpressureA = 0;
	extern double vesselpressureA;
	vesselpressureA = 0;
		extern double hompressureA;
	hompressureA = 0;
	VFWater = RetrieveDoubleParameter( ParameterFile, "VFWater" );
	VFCell = RetrieveDoubleParameter( ParameterFile, "VFCell" );
	ChMobility = RetrieveDoubleParameter( ParameterFile, "ChMobility" );
	ChLength = RetrieveDoubleParameter( ParameterFile, "ChLength" );
	MitosisRate = RetrieveDoubleParameter(ParameterFile, "MitosisRate");
 
 extern bool MinimalSaveMode;
 MinimalSaveMode = (bool) RetrieveDoubleParameter( ParameterFile, "MinimalSaveMode" );
 
 extern char* StringBuffer;
 extern int StringBufferLength;
 extern char* LogFile;

 RetrieveStringParameter( ParameterFile , "SaveDirectory" ,SaveDirectory );
 ResumeFromPreviousData = (bool) RetrieveDoubleParameter( ParameterFile, "ResumeFromPreviousData" );
 InitializeSimulationFromPictures = 
     (bool) RetrieveDoubleParameter( ParameterFile, "InitializeSimulationFromPictures" );
 ScratchSpaceMode = (bool) RetrieveDoubleParameter( ParameterFile, "ScratchSpaceMode" );
 
 IsFirstOxygenUpdate = true;
 int i,j;
 // initialize the 
 
 CreateFilebase( PhiFile, "Phi_resume.dat" );
 InsertDirectory( PhiFile, SaveDirectory );
 CreateFilebase( Xfile, "X_resume.dat" );
 InsertDirectory( Xfile, SaveDirectory );
 CreateFilebase( Yfile, "Y_resume.dat" );
 InsertDirectory( Yfile, SaveDirectory );
 CreateFilebase( PhiNecroticFile , "PhiN_resume.dat" );
 InsertDirectory( PhiNecroticFile, SaveDirectory );
 CreateFilebase( PhiQuiescentFile , "PhiH_resume.dat" );
 InsertDirectory( PhiQuiescentFile, SaveDirectory );
 CreateFilebase( OxygenFile , "Oxygen_resume.dat" );
 InsertDirectory( OxygenFile, SaveDirectory );
 CreateFilebase( OxygenDiffusivityFile , "OxygenDiffusivity_resume.dat" );
 InsertDirectory( OxygenDiffusivityFile, SaveDirectory );
 CreateFilebase( PressureFile , "Pressure_resume.dat" );
 InsertDirectory( PressureFile, SaveDirectory );
 CreateFilebase( MobilityFile , "Mobility_resume.dat" );
 InsertDirectory( MobilityFile, SaveDirectory );

 CreateFilebase( BloodVesselDensityFile , "BloodVesselDensity_resume.dat" );
 InsertDirectory( BloodVesselDensityFile, SaveDirectory );
 
 CreateFilebase( ECM_File , "ECM_resume.dat" );
 InsertDirectory( ECM_File, SaveDirectory );
 CreateFilebase( MDEfile , "MDE_resume.dat" ); 
 InsertDirectory( MDEfile, SaveDirectory );
 CreateFilebase( TAFfile , "TAF_resume.dat" ); 
 InsertDirectory( TAFfile, SaveDirectory ); 
 
 CreateFilebase( StructureFile , "Structure_resume.dat" );
 InsertDirectory( StructureFile , SaveDirectory );

 CreateFilebase( TumorParameterFile , "Parameters_resume.cfg" );
 InsertDirectory( TumorParameterFile , SaveDirectory );
 
 CreateFilebase( BloodVesselIndicatorXFile , "BloodVesselIndicatorX_resume.dat" );
 InsertDirectory( BloodVesselIndicatorXFile , SaveDirectory );
 CreateFilebase( BloodVesselIndicatorYFile , "BloodVesselIndicatorY_resume.dat" );
 InsertDirectory( BloodVesselIndicatorYFile , SaveDirectory );
 CreateFilebase( BloodVesselIndicatorFile , "BloodVesselIndicator_resume.dat" ); 
 InsertDirectory( BloodVesselIndicatorFile , SaveDirectory );
 CreateFilebase( HaematocritLevelFile , "Haematocrit_resume.dat" );
 InsertDirectory( HaematocritLevelFile, SaveDirectory );
  CreateFilebase( BloodVesselRadiusXdirFile , "BloodVesselRadiusXdir_resume.dat" );
 InsertDirectory( BloodVesselRadiusXdirFile, SaveDirectory );
  CreateFilebase( BloodVesselRadiusYdirFile , "BloodVesselRadiusYdir_resume.dat" );
 InsertDirectory( BloodVesselRadiusYdirFile, SaveDirectory );
  CreateFilebase( BloodVesselRadiusFile , "BloodVesselRadius_resume.dat" );
 InsertDirectory( BloodVesselRadiusFile, SaveDirectory ); 
 CreateFilebase(SproutTipsFile , "SproutTips_resume.dat" );
 InsertDirectory( SproutTipsFile, SaveDirectory );
 CreateFilebase( SproutAgeFile , "SproutAge_resume.dat" ); 
 InsertDirectory( SproutAgeFile, SaveDirectory );
 
 CreateFilebase( ConductXFile , "ConductX_resume.dat" );
 InsertDirectory( ConductXFile, SaveDirectory );
 CreateFilebase( ConductYFile , "ConductY_resume.dat" );
 InsertDirectory( ConductYFile, SaveDirectory );

 CreateFilebase( VesselAgeFile , "VesselAge_resume.dat" ); 
 InsertDirectory( VesselAgeFile, SaveDirectory ); 


 
 Xleft = RetrieveDoubleParameter( ParameterFile, "Xleft" );
 Xright = RetrieveDoubleParameter( ParameterFile, "Xright" );
 Xnodes = (int) RetrieveDoubleParameter( ParameterFile, "Xnodes" );

 Ybottom = RetrieveDoubleParameter( ParameterFile, "Ybottom" );
 Ytop = RetrieveDoubleParameter( ParameterFile, "Ytop" );
 Ynodes = (int) RetrieveDoubleParameter( ParameterFile, "Ynodes" );
 
 SimulationStartTime = RetrieveDoubleParameter( ParameterFile, "SimulationStartTime" );
 SimulationEndTime = RetrieveDoubleParameter( ParameterFile, "SimulationEndTime" );
 OutputInterval = RetrieveDoubleParameter( ParameterFile, "OutputInterval" );

 TumorOxygenUptakeRate = RetrieveDoubleParameter( ParameterFile, "TumorOxygenUptakeRate" );
 QuiescentOxygenUptakeRate = RetrieveDoubleParameter( ParameterFile, "QuiescentOxygenUptakeRate" );
 NecroticOxygenUptakeRate = RetrieveDoubleParameter( ParameterFile, "NecroticOxygenUptakeRate" );
 ECMoxygenUptakeRate = RetrieveDoubleParameter( ParameterFile, "ECMoxygenUptakeRate" );
 ExteriorOxygenDiffusivity = RetrieveDoubleParameter( ParameterFile, "ExteriorOxygenDiffusivity" );
 AdvectionConstant = RetrieveDoubleParameter( ParameterFile, "AdvectionConstant" );
 bulk_blood_tissue_transfer_rate = RetrieveDoubleParameter( ParameterFile, "bulk_blood_tissue_transfer_rate" );
 bulk_blood_transfer_cutoff_distance = 
    RetrieveDoubleParameter( ParameterFile, "bulk_blood_transfer_cutoff_distance" );
 TopBloodVessel = (bool) RetrieveDoubleParameter( ParameterFile, "TopBloodVessel" );
 BottomBloodVessel = (bool) RetrieveDoubleParameter( ParameterFile, "BottomBloodVessel" );
    
 // MaxOxygenDeliveryRate = RetrieveDoubleParameter( ParameterFile, "MaxOxygenDeliveryRate" );

 QuiescenceOxygenLevel = RetrieveDoubleParameter( ParameterFile, "QuiescenceOxygenLevel" );
 NecrosisOxygenLevel = RetrieveDoubleParameter( ParameterFile, "NecrosisOxygenLevel" );
 ExteriorMobility = RetrieveDoubleParameter( ParameterFile, "ExteriorMobility" );
 
 ApoptosisRate = RetrieveDoubleParameter( ParameterFile, "ApoptosisRate" );
 TumorAggressivenessG = RetrieveDoubleParameter( ParameterFile, "TumorAggressivenessG" );
 TumorNecrosisRateGN = RetrieveDoubleParameter( ParameterFile, "TumorNecrosisRateGN" );
 DRG_Effect = RetrieveDoubleParameter( ParameterFile, "DRG_Effect" );

 
 ECM_DegradationConstant = RetrieveDoubleParameter( ParameterFile, "ECM_DegradationConstant" );
 ECM_ProductionConstant = RetrieveDoubleParameter( ParameterFile, "ECM_ProductionConstant" );
 
 TAFdiffusionConstant = RetrieveDoubleParameter( ParameterFile, "TAFdiffusionConstant" );
 TAFprodctionConstant = RetrieveDoubleParameter( ParameterFile, "TAFprodctionConstant" );
 TAFdecayConstant = RetrieveDoubleParameter( ParameterFile, "TAFdecayConstant" );
 
 EndothelialTAFuptakeRate = RetrieveDoubleParameter( ParameterFile, "EndothelialTAFuptakeRate" );
 
 MDEdiffusionConstant = RetrieveDoubleParameter( ParameterFile, "MDEdiffusionConstant" );
 MaxMDEdiffusionConstant = RetrieveDoubleParameter( ParameterFile, "MaxMDEdiffusionConstant" );
 if( MaxMDEdiffusionConstant < MDEdiffusionConstant )
 { MaxMDEdiffusionConstant = 1.0; }
 
 // create the MDE polynomial;
 CreateMDEdiffusivityPolynomial();
 
 
 MDEprodctionConstant = RetrieveDoubleParameter( ParameterFile, "MDEprodctionConstant" );
 MDEdecayConstant = RetrieveDoubleParameter( ParameterFile, "MDEdecayConstant" );
 MDE_DegradationConstant = RetrieveDoubleParameter( ParameterFile, "MDE_DegradationConstant" );
 
 DRGprodctionConstant = RetrieveDoubleParameter( ParameterFile, "DRGprodctionConstant" );
 DRGDecayVessel = RetrieveDoubleParameter( ParameterFile, "DRGDecayVessel" );
 DRGdecayConstant = RetrieveDoubleParameter( ParameterFile, "DRGdecayConstant" );
 DRG_DegradationConstant = RetrieveDoubleParameter( ParameterFile, "DRG_DegradationConstant" );
 D_DRG = RetrieveDoubleParameter( ParameterFile, "D_DRG" );

 
 TumorProliferationConstant = RetrieveDoubleParameter( ParameterFile, "TumorProliferationConstant" );
 TumorSurfaceTensionConstant = RetrieveDoubleParameter( ParameterFile, "TumorSurfaceTensionConstant" );
 Custom_Cutoff_Max = RetrieveDoubleParameter( ParameterFile , "Custom_Cutoff_Max" ); 
 
 collapseConstant = RetrieveDoubleParameter( ParameterFile , "collapseConstant" ); 
 extern double MinimumMobility;
 MinimumMobility = RetrieveDoubleParameter(ParameterFile,"MinimumMobility");
 if( MinimumMobility < 1e-8 )
 { MinimumMobility = 1e-8; }
 
 extern double MaximumMobility;
 MaximumMobility = RetrieveDoubleParameter(ParameterFile,"MaximumMobility");
 
 extern double MobilityECM_ParameterA;
 MobilityECM_ParameterA = 1.0/MinimumMobility - 1.0/MaximumMobility;
 
 extern double MobilityECM_ParameterB;
 MobilityECM_ParameterB = 1.0/MaximumMobility;
 
 extern bool FunctionalMobility;
 FunctionalMobility = (bool) RetrieveDoubleParameter(ParameterFile,"FunctionalMobility");
 
 extern double FunctionSmoothingLevel;
 FunctionSmoothingLevel = RetrieveDoubleParameter(ParameterFile,"FunctionSmoothingLevel");
 
 
 if( ResumeFromPreviousData )
 {
  InitializeSimulationFromPictures = false;
 
  StartStringBuffer( "Resuming from saved data in directory " );
  AddToStringBuffer( SaveDirectory );
  AddToStringBuffer( " ... \n" );
  cout << StringBuffer << endl;
  MakeLogEntry( LogFile, StringBuffer );
  
  X.ReadFromFile( Xfile );
  Y.ReadFromFile( Yfile );

  Xnodes = X.TellCols();
  Xleft = *X(0);
  Xright = *X(Xnodes-1);

  Ynodes = Y.TellCols();
  Ybottom = *Y(0);
  Ytop = *Y(Ynodes-1);

  Phi.ReadFromFile( PhiFile );
  PhiNecrotic.ReadFromFile( PhiNecroticFile );
  PhiQuiescent.ReadFromFile( PhiQuiescentFile );
  if( PhiQuiescent.TellRows() < Xnodes || PhiQuiescent.TellCols() < Ynodes )
  { PhiQuiescent = ones(Xnodes,Ynodes); 
  cout<<"PhiHlost"<<endl;} 
 
  MDE.ReadFromFile( MDEfile );
  if( MDE.TellRows() < Xnodes || MDE.TellCols() < Ynodes )
  { MDE = zeros(Xnodes,Ynodes); 
  cout<<"MDE LOST"<<endl;}  
  
  TAF.ReadFromFile( TAFfile );
  if( TAF.TellRows() < Xnodes || TAF.TellCols() < Ynodes )
  { TAF = zeros(Xnodes,Ynodes); }  
  
  ECM.ReadFromFile( ECM_File );
  if( ECM.TellRows() < Xnodes || ECM.TellCols() < Ynodes )
  { ECM = ones(Xnodes,Ynodes); }
  
  
  char Original_ECM_File [1024];
  sprintf( Original_ECM_File , "%s/OriginalECM_resume.dat" , SaveDirectory );
  Original_ECM.ReadFromFile( Original_ECM_File );
  if( Original_ECM.TellRows() < Xnodes || Original_ECM.TellCols() < Ynodes )
  { Original_ECM = ECM; }
  
  Oxygen.ReadFromFile( OxygenFile );
  if( Oxygen.TellRows() < Xnodes || Oxygen.TellCols() < Ynodes )
  { Oxygen = ones(Xnodes,Ynodes); }

  Pressure.ReadFromFile( PressureFile );
  if( Pressure.TellRows() < Xnodes || Pressure.TellCols() < Ynodes )
  { Pressure = zeros(Xnodes,Ynodes); }
  
  Mobility.ReadFromFile( MobilityFile );
  if( Mobility.TellRows() < Xnodes || Mobility.TellCols() < Ynodes )
  { Mobility = ones(Xnodes,Ynodes); }
  
  OxygenDiffusivity.ReadFromFile( OxygenDiffusivityFile );
  if( OxygenDiffusivity.TellRows() < Xnodes || OxygenDiffusivity.TellCols() < Ynodes )
  { OxygenDiffusivity = ones(Xnodes,Ynodes); }
  
  BloodVesselIndicatorX.ReadFromFile( BloodVesselIndicatorXFile );
  if( BloodVesselIndicatorX.TellRows() < Xnodes || BloodVesselIndicatorX.TellCols() < Ynodes )
  { BloodVesselIndicatorX = zeros(Xnodes,Ynodes); 
   cout<<"BVIX information lost"<<endl;} 
   
  BloodVesselIndicatorY.ReadFromFile( BloodVesselIndicatorYFile );
  if( BloodVesselIndicatorY.TellRows() < Xnodes || BloodVesselIndicatorY.TellCols() < Ynodes )
  { BloodVesselIndicatorY = zeros(Xnodes,Ynodes); 
   cout<<"BVIY information lost"<<endl;} 
     

  BloodVesselIndicator.ReadFromFile( BloodVesselIndicatorFile );
  if( BloodVesselIndicator.TellRows() < Xnodes || BloodVesselIndicator.TellCols() < Ynodes )
  { BloodVesselIndicator = zeros(Xnodes,Ynodes); } 
 
   SproutTips.ReadFromFile( SproutTipsFile );
  //cout<<"say something"<<endl;
  if( SproutTips.TellRows() < Xnodes || SproutTips.TellCols() < Ynodes )
  { SproutTips = zeros(Xnodes,Ynodes); 
  cout<<"26say something about 0"<<endl;} 
  
  SproutAge.ReadFromFile( SproutAgeFile );
  if( SproutAge.TellRows() < Xnodes || SproutAge.TellCols() < Ynodes )
  { SproutAge = zeros(Xnodes,Ynodes); 
  cout<<"sproutAge information lost, set to 0"<<endl;}
  
  VesselAge.ReadFromFile( VesselAgeFile );
  if( VesselAge.TellRows() < Xnodes || VesselAge.TellCols() < Ynodes )
  { VesselAge = zeros(Xnodes,Ynodes);
  cout<<"VesselAge information lost, set to 0"<<endl; }


  Structure.ReadFromFile( StructureFile );
  if( Structure.TellRows() < Xnodes || Structure.TellCols() < Ynodes )
  { Structure = zeros(Xnodes,Ynodes); } 
  
  BloodVesselDensity.ReadFromFile( BloodVesselDensityFile );
  if( BloodVesselDensity.TellRows() < Xnodes || BloodVesselDensity.TellCols() < Ynodes )
  { BloodVesselDensity = ones(Xnodes,Ynodes); } 
  
  HaematocritLevel.ReadFromFile( HaematocritLevelFile );
  if( HaematocritLevel.TellRows() < Xnodes ||HaematocritLevel.TellCols() < Ynodes )
  { HaematocritLevel = zeros(Xnodes,Ynodes); 
  cout<<"Haematocritlevel information lost, set to 0"<<endl;} 
  
  BloodVesselRadiusXdir.ReadFromFile( BloodVesselRadiusXdirFile );
  if( BloodVesselRadiusXdir.TellRows() < Xnodes || BloodVesselRadiusXdir.TellCols() < Ynodes )
  { 
  
  for(int i=0;i<Xnodes;i++)
  {
	for(int j=0;j<Ynodes;j++)
	{
		*BloodVesselRadiusXdir(i,j) = 1e-6;
	}
  
  }
  cout<<"RADIUS information lost, set to uniform 1e-6"<<endl; } 
  BloodVesselRadiusYdir.ReadFromFile( BloodVesselRadiusYdirFile );
  if(BloodVesselRadiusYdir.TellRows() < Xnodes || BloodVesselRadiusYdir.TellCols() < Ynodes )
  { 
  for(int i=0;i<Xnodes;i++)
  {
	for(int j=0;j<Ynodes;j++)
	{
		*BloodVesselRadiusYdir(i,j) = 1e-6;
	}
  
  }
  
  cout<<"RADIUS information lost, set to uniform 1e-6"<<endl;
   } 
  BloodVesselRadius.ReadFromFile( BloodVesselRadiusFile );
  if(BloodVesselRadius.TellRows() < Xnodes || BloodVesselRadius.TellCols() < Ynodes )
  { BloodVesselRadius = zeros(Xnodes,Ynodes);
 // cout<<"25say something about 0"<<endl; 
   } 
  
  
   ConductX.ReadFromFile( ConductXFile );
  if( ConductX.TellRows() < Xnodes || ConductX.TellCols() < Ynodes )
  { ConductX = zeros(Xnodes,Ynodes);
  //cout<<"21say something about 0"<<endl;
   } 
  ConductY.ReadFromFile( ConductYFile );
  if(ConductY.TellRows() < Xnodes || ConductY.TellCols() < Ynodes )
  { ConductY = zeros(Xnodes,Ynodes); }
	 
	
  
  
  
  SimulationStartTime = RetrieveDoubleParameter( TumorParameterFile, "time" );
  cout<<SimulationStartTime<<endl;
  TumorOxygenUptakeRate = RetrieveDoubleParameter( TumorParameterFile, "TumorOxygenUptakeRate" );
  QuiescentOxygenUptakeRate = RetrieveDoubleParameter( TumorParameterFile, "QuiescentOxygenUptakeRate" );
  NecroticOxygenUptakeRate = RetrieveDoubleParameter( TumorParameterFile, "NecroticOxygenUptakeRate" );
  ECMoxygenUptakeRate = RetrieveDoubleParameter( TumorParameterFile, "ECMoxygenUptakeRate" );
//  MaxOxygenDeliveryRate = RetrieveDoubleParameter( TumorParameterFile, "MaxOxygenDeliveryRate" );

  QuiescenceOxygenLevel = RetrieveDoubleParameter( TumorParameterFile, "QuiescenceOxygenLevel" );
  NecrosisOxygenLevel = RetrieveDoubleParameter( TumorParameterFile, "NecrosisOxygenLevel" );
 
  ApoptosisRate = RetrieveDoubleParameter( TumorParameterFile, "ApoptosisRate" );
  TumorAggressivenessG = RetrieveDoubleParameter( TumorParameterFile, "TumorAggressivenessG" );
  TumorNecrosisRateGN = RetrieveDoubleParameter( TumorParameterFile, "TumorNecrosisRateGN" );
  
 }
 else
 {
 
  
  X = linspace( Xleft , Xright, Xnodes );
  Y = linspace( Ybottom, Ytop, Ynodes );



  Phi = zeros( Xnodes, Ynodes );
 // getchar();
  PhiNecrotic = ones( Xnodes,Ynodes);
  PhiQuiescent = ones( Xnodes, Ynodes );

  MDE = zeros( Xnodes, Ynodes );
  Macrof = zeros(Xnodes, Ynodes);
  TracerT = zeros(Xnodes, Ynodes);
  TAF = zeros( Xnodes, Ynodes );
  ECM = ones( Xnodes, Ynodes );
  TracerD = zeros(Xnodes, Ynodes);
  Debuger = zeros(Xnodes, Ynodes);
  Oxygen = ones( Xnodes, Ynodes );
  OxygenDiffusivity = ones( Xnodes , Ynodes );
  
  Pressure = zeros( Xnodes, Ynodes );
  Mobility = ones( Xnodes, Ynodes );
  SrtX = zeros( Xnodes, Ynodes );
  SrtY = zeros( Xnodes, Ynodes );
	 PressureGradientX = zeros( Xnodes, Ynodes );
	 PressureGradientY = zeros( Xnodes, Ynodes );
	 Swss = zeros( Xnodes, Ynodes);
	 Sm = zeros( Xnodes, Ynodes);
	 Sp = zeros( Xnodes, Ynodes);
	 Ss = zeros( Xnodes, Ynodes);
	 Svar = zeros( Xnodes, Ynodes);
	 
	 SwssX = zeros( Xnodes, Ynodes);
	 SmX = zeros( Xnodes, Ynodes);
	 SpX = zeros( Xnodes, Ynodes);
	 SsX = zeros( Xnodes, Ynodes);
	 SvarX = zeros( Xnodes, Ynodes);
	 SwssY= zeros( Xnodes, Ynodes);
	 SmY = zeros( Xnodes, Ynodes);
	 SpY = zeros( Xnodes, Ynodes);
	 SsY = zeros( Xnodes, Ynodes);
	 SvarY = zeros( Xnodes, Ynodes);
  TotalPressure = zeros(Xnodes,Ynodes);

  BloodVesselIndicator = zeros(Xnodes,Ynodes);
  BloodVesselIndicatorX = zeros(Xnodes,Ynodes);
  BloodVesselIndicatorY = zeros(Xnodes,Ynodes);
  VEGFR = zeros(Xnodes, Ynodes);
  Notch = zeros(Xnodes, Ynodes);
  Delta4 = zeros(Xnodes, Ynodes);
  VesselCells = zeros(Xnodes, Ynodes);
  Shuffle = zeros(Xnodes, Ynodes);
  HomPressure = zeros(Xnodes,Ynodes);
  Tracer = zeros(Xnodes, Ynodes);
  TracerS = zeros(Xnodes, Ynodes);
  HaematocritLevel = zeros(Xnodes, Ynodes);
  BloodVesselDensity = ones(Xnodes,Ynodes);
  SproutTips = zeros( Xnodes, Ynodes );
  Structure = zeros( Xnodes, Ynodes);
  CollapseClock = zeros( Xnodes, Ynodes);
  double CenterX = 0.5*(Xleft + Xright);
  double CenterY = 0.5*(Ybottom + Ytop );
  double Radius = MIN( Ytop - CenterY, Xright - CenterX )/2.0;
    
  for( i=0 ; i < Xnodes ; i++ )
  {
   for( j=0 ; j < Ynodes ; j++ )
   {
    *Phi(i,j) = sqrt( square(*X(i)-CenterX)+square(*Y(j)-CenterY) ) - Radius;
	*HomPressure(i,j) = 2000;
   }
  }
  
 }
 
 dx = *X(1) - *X(0);
 dy = *Y(1) - *Y(0);
 dt = 0.00;
 t = SimulationStartTime;
 
 CreateO2uptakeFunction();

  
 if( InitializeSimulationFromPictures )
 {
  RetrieveStringParameter( ParameterFile , "TumorPictureFile"  ,TumorPictureFile );
  RetrieveStringParameter( ParameterFile , "TissuePictureFile" ,TissuePictureFile );
  
  RetrieveStringParameter( ParameterFile , "TissueTypesFile", TissueTypesFile );
 
  StartStringBuffer( "Note: Initializing tumor from picture " );
  AddToStringBuffer( TumorPictureFile );
  AddToStringBuffer( " ... \n" );
  cout << StringBuffer;
  MakeLogEntry( LogFile, StringBuffer );
  
  InitializeTumorFromPicture( TumorPictureFile , TissuePictureFile );
  /***
  StartStringBuffer( "Note: Initializing tissue from picture " );
  AddToStringBuffer( TissuePictureFile );
  AddToStringBuffer( " ... \n" );
  cout << StringBuffer;
  MakeLogEntry( LogFile, StringBuffer );
  
  //InitializeTissueFromPicture( TissuePictureFile, TissueTypesFile );

  cout << "\t" << "oxygen" << endl;
  UpdateOxygen();
  cout << "\t" << "taf" << endl;
  UpdateGlobalTAF();
  cout << "\t" << "microenvironment" << endl;
  UpdateMicroenvironment();
  cout << "\t" << "pressure" << endl;
  UpdatePressure();
****/
 }
 
 int k = 0;
 while( k*OutputInterval <= t + _EPS )
 { k++; }
 NextOutputTime = k*OutputInterval;
 if( fabs( t - (k-1)*OutputInterval ) <= _EPS )
 { NextOutputTime = (k-1)*OutputInterval; }
 
 //MaxMobility = Mobility.InfinityNorm();
 
 //SaveTumorData( t ); 
 //SaveResumeTumorData();
 cout<<"before"<<endl;
 //getchar();

 //QuickPlot();
 cout<<"after"<<endl;
 //getchar();
 
 return true;
}

bool InitializeTumorFromPicture( const char* FileName , const char* StructureFileName )
{
 extern double Xleft;
 extern double Xright;
 extern double Ybottom;
 extern double Ytop;
 extern int Xnodes;
 extern int Ynodes;
 extern bool TumorBegin;
 extern Matrix Phi;
 extern Matrix PhiNecrotic;
 extern Matrix PhiQuiescent;
 extern Matrix X;
 extern Matrix Y;
 extern double dx;
 extern double dy;
 extern double dt;
 
 extern double _EPS;
 
 extern char* ParameterFile;
 extern char SaveDirectory [1025];
 extern char* StringBuffer;
 extern char* LogFile;
 
 extern Matrix MDE;
 extern Matrix TAF;
 extern Matrix ECM;
 extern Matrix Structure;
 
 extern Matrix Oxygen;
 extern Matrix Pressure;
 extern Matrix Mobility;
 extern Matrix OxygenDiffusivity;
 
 extern Matrix BloodVesselIndicator;
 extern Matrix BloodVesselRadius;
 extern Matrix Flow;
 extern Matrix SproutTips;
 
 int i,j;
 
// BMP Gray;
// Gray.ReadFromFile( FileName );

// BMP ReferenceImage;
// ReferenceImage.ReadFromFile( StructureFileName );
/**
 if( ReferenceImage.TellWidth() != Gray.TellWidth() ||
     ReferenceImage.TellHeight() != Gray.TellHeight() )
 { 
  ReferenceImage.SetSize( Gray.TellWidth() , Gray.TellHeight() );
  
  for( int jj=0; jj < Gray.TellHeight() ;jj++ )
  {
   for( int ii=0 ; ii < Gray.TellWidth() ; ii++ )
   {
    *ReferenceImage(i,j) = *Gray(0,0);
   }
  }
 
 }
 **/
 
// ImageDifference( ReferenceImage, Gray );
 
 extern double MicronsPerPixel;
 extern double MicronsPerDeltaX;
  
 MicronsPerPixel = RetrieveDoubleParameter(ParameterFile,"MicronsPerPixel");
 MicronsPerDeltaX = RetrieveDoubleParameter(ParameterFile,"MicronsPerDeltaX");
 
 extern double DiffusionLengthScaleInMicrons;
 DiffusionLengthScaleInMicrons = RetrieveDoubleParameter(ParameterFile,"DiffusionLengthScaleInMicrons");
 if( DiffusionLengthScaleInMicrons < 1 ) 
 { DiffusionLengthScaleInMicrons = 200.0; } // default value of 200 microns
 
 // first, get the nondimensional domain width and height. 
 // Then, determine the number of x and y nodes required.
 // The, discretize X and Y.
 
// double domain_width = ( Gray.TellWidth() - 1.0 ) * MicronsPerPixel / DiffusionLengthScaleInMicrons;
 Xleft = -10;//-10;
 Xright = 10;//10;
 dx = 0.1;//0.10; // 1/10 of length scale is appropriate
 Xnodes = (int) ceil( (Xright-Xleft)/dx + 1.0 );
// Xnodes = (int) ImprovedRound( (Xright-Xleft)/dx + 1.0 );
 X = linspace( Xleft , Xright, Xnodes );
 //cout<<*X(1)<<endl;
 
// double domain_height = ( Gray.TellHeight() - 1.0 ) * MicronsPerPixel / DiffusionLengthScaleInMicrons;
 Ybottom = -10;//-10;
 Ytop = 10;//10;
 dy = dx = 0.1;//0.10; // 1/10 of length scale is appropriate
 Ynodes = (int) ceil( (Ytop-Ybottom)/dy + 1.0 );
// Ynodes = (int) ImprovedRound( (Ytop-Ybottom)/dy + 1.0 );
 Y = linspace( Ybottom , Ytop , Ynodes );
//cout<<*Y(1)<<endl;
 // rescale the image. 
 
// Rescale( Gray, 'W', Xnodes );
 
// cout << "xn,yn: " << Xnodes << " x " << Ynodes << endl
 //     << "image: " << Gray.TellWidth() << " x " << Gray.TellHeight() << endl;
/** 
 Phi = BMPtoMatrix( Gray,-1,1 );
 if( Phi.TellRows() > Xnodes || Phi.TellCols() > Ynodes )
 {
  Phi = Phi.submatrix( 0,Xnodes-1,0,Ynodes-1 );
 }
 
// if( Phi.TellRows() != Xnodes || Phi.TellCols() != Ynodes )
// {
  StartStringBuffer( "Warning: tumor image wrong size! Creating circular tumor" );
  AddToStringBuffer( SaveDirectory );
  AddToStringBuffer( " ... \n" );
  cout << StringBuffer << endl;
  MakeLogEntry( LogFile, StringBuffer );
 **/
  Phi = ones(Xnodes,Ynodes);
  
  for( i=0 ; i < Xnodes ; i++ )
  {
   for( int j=0; j < Ynodes ; j++ )
   {
   // *Phi(i,j) = sqrt( Square(i - (Xnodes-1)/2.0) + Square(j-(Ynodes-1)/2.0) ) - 0.045*Xnodes;
    *Phi(i,j) = sqrt( Square(i - (Xnodes-1)/2.0) + Square(j-(Ynodes-1)/2.0) ) - 0.025*Xnodes;
   }
  }
  
// }


 Phi = GaussianFilter( X,Y,Phi,2*dx ); // NonExtended
//  cout<<"after filter"<<endl;
// getchar();
 Phi = DiscreteReinit( Phi );
 Phi = Reinit( Phi );
 

 GenerateGeometricGlobals( X,Y,Phi );
  ECM = ones(Xnodes,Ynodes);

 MDE = zeros(Xnodes,Ynodes);
 TAF = zeros(Xnodes,Ynodes);
 BloodVesselIndicator = zeros(Xnodes,Ynodes);
 BloodVesselRadius = zeros(Xnodes,Ynodes);
 Flow = zeros(Xnodes,Ynodes);
 SproutTips = zeros(Xnodes,Ynodes);

 Oxygen = ones(Xnodes,Ynodes);
 OxygenDiffusivity = ones(Xnodes,Ynodes);
 
 bool ReinitQ = false;
 bool ReinitN = false;
 /*
 for( int j=0; j < Ynodes ; j++ )
 {
  for( int i=0; i < Xnodes ; i++ )
  {
   if( *PhiQuiescent(i,j) <= _EPS ){ ReinitQ = true; cout<< *PhiQuiescent(i,j);}
   if( *PhiNecrotic(i,j) <= _EPS ){ ReinitN = true; }
  }
 } 
 */
 
 if( ReinitQ )
 {
  PhiQuiescent = Reinit( PhiQuiescent );
  PhiQuiescent = Reinit( PhiQuiescent );
  PhiQuiescent = Reinit( PhiQuiescent );

  PhiQuiescent = Reinit( PhiQuiescent );
  PhiQuiescent = Reinit( PhiQuiescent );
  PhiQuiescent = Reinit( PhiQuiescent );
 }
 if( ReinitN )
 {
  PhiNecrotic = Reinit( PhiNecrotic );
  PhiNecrotic = Reinit( PhiNecrotic );
  PhiNecrotic = Reinit( PhiNecrotic ); 

  PhiNecrotic = Reinit( PhiNecrotic );
  PhiNecrotic = Reinit( PhiNecrotic );
  PhiNecrotic = Reinit( PhiNecrotic );
 }

 Pressure = zeros(Xnodes,Ynodes);
 Mobility = ones(Xnodes,Ynodes);
 ECM = ones(Xnodes,Ynodes);
  
 TAF = zeros(Xnodes,Ynodes);
 MDE = zeros(Xnodes,Ynodes);
 
 Structure = zeros(Xnodes,Ynodes);
 dt = 0;
 TumorBegin = RetrieveDoubleParameter(ParameterFile,"TumorBegin");
 return true;
}

bool InitializeBloodVesselsFromPicture( const char* FileName )
{
 extern char* ParameterFile;
 extern char SaveDirectory [1025];
 extern char* StringBuffer;
 extern char* LogFile;

 extern int Xnodes;
 extern int Ynodes;
 
 extern Matrix BloodVesselIndicator;
 extern Matrix BloodVesselRadius;
 extern Matrix Flow;

 int i,j;

 BMP Gray;
 Gray.ReadFromFile( FileName );
 
 extern double MicronsPerPixel;
 extern double MicronsPerDeltaX;
  
 MicronsPerPixel = RetrieveDoubleParameter(ParameterFile,"MicronsPerPixel");
 MicronsPerDeltaX = RetrieveDoubleParameter(ParameterFile,"MicronsPerDeltaX");

 extern double DiffusionLengthScaleInMicrons;
 DiffusionLengthScaleInMicrons = RetrieveDoubleParameter(ParameterFile,"DiffusionLengthScaleInMicrons");
 if( DiffusionLengthScaleInMicrons < 1 ) 
 { DiffusionLengthScaleInMicrons = 200.0; } // default value of 200 microns

 Rescale( Gray, 'W', Xnodes );
 
 BloodVesselIndicator = BMPtoMatrix( Gray , 1 , 0 );
 if( BloodVesselIndicator.TellRows() > Xnodes ||BloodVesselIndicator.TellCols() > Ynodes )
 {
  BloodVesselIndicator = BloodVesselIndicator.submatrix( 0,Xnodes-1,0,Ynodes-1 );
 } 
 
 if( BloodVesselIndicator.TellRows() != Xnodes ||
     BloodVesselIndicator.TellCols() != Ynodes )
 {
  StartStringBuffer( "Warning: blood vessel image wrong size! Using default bv == 0" );
  AddToStringBuffer( SaveDirectory );
  AddToStringBuffer( " ... \n" );
  cout << StringBuffer << endl;
  MakeLogEntry( LogFile, StringBuffer );
 
  BloodVesselIndicator = zeros(Xnodes,Ynodes);
  return false;
 }   
 
 for( int i=0 ; i < Xnodes ; i++ )
 {
  for( int j=0; j < Ynodes ; j++ )
  {
   double bvi = *BloodVesselIndicator(i,j);
   if( bvi > 0.01 && bvi < 0.5 )
   { *BloodVesselIndicator(i,j) = 0.0; }
   if( bvi > 0.49 && bvi < 0.99 )
   { *BloodVesselIndicator(i,j) = 1.0; }
  }  
 } 
 
 BloodVesselRadius = zeros(Xnodes,Ynodes);
 Flow = zeros(Xnodes,Ynodes);
 
 double TempRadius = RetrieveDoubleParameter( ParameterFile , "Rmin" );

 for( i=0 ; i < Xnodes ; i++ )
 {
  for( j=0 ; j < Ynodes ; j++ )
  {
   if( *BloodVesselIndicator(i,j) > 0.25 )
   { *BloodVesselRadius(i,j) = TempRadius; }
  }
 } 

 return true;
}

bool InitializeTissueFromPicture( const char* ImageFile , const char* TissueCharacteristicsFile )
{
 extern Matrix Structure;
 extern Matrix ECM;
 extern Matrix BloodVesselDensity;
 extern Matrix Phi;
 extern Matrix Mobility;
 extern int Xnodes;
 extern int Ynodes;

 bool return_value = InitializeTissueTypesFromXML( TissueCharacteristicsFile );
 DisplayAllTissueTypes();
 
 return_value &= ImportTissues( ImageFile );

 //Structure.WriteToTextFile("S.dat");
 //ECM.WriteToTextFile("ECM.dat");
 //BloodVesselDensity.WriteToTextFile("BVD.dat");
 //Phi.WriteToTextFile( "Phi.dat");
 cout<<ECM.TellRows()<<endl;
  cout<<"Xnodes"<<Xnodes<<"Ynodes"<<Ynodes<<"ECM 201 201"<<*ECM(201,201)<<endl;
 //getchar(); 
 //UpdateMobility();
 
 //Mobility.WriteToTextFile( "Mu.dat");
 
 return return_value;
}
